DALS004-统计推断(Inference)02-CLT与T-test

发表于 | 分类于
| 字数统计: 9,886 | 阅读时长 ≈ 39

前言

这篇笔记是《Data Analysis for the Life Sciences》的第2章统计推断(Inference)的第3部分。这一部分的主要内容涉及中心极限理论，t检验和功效计算。

中心极限定理(Central Limit Theorem)与t分布

中心极限定理(CLT)

当样本的数目非常大时，一个随机样本的均值 $\bar{Y}$ 就会服从正态分布，这个正态分布的均值是总体均值( $\mu_{Y}$ )，标准差是总体的标准差 $\sigma_{Y}$ 除以样本的数目 $N$ 的平方根。

如果我们将一个随机变量减去一个常数，那么这个新生成的随机变量就由该常数转换而来的（原文是：Please note that if we subtract a constant from a random variable, the mean of the new random variable shifts by that constant。我个人的理解就是，这个新的变量就是与原来的变量相差这个常数），从数学角度来说，如果$X$是一个服从均值为 $\mu$ 的分布，$a$是一个常数，那么 $X-a$ 的均值的分布的均值就是$\mu-a$。这个也适应用乘法与标准差(SD)。

如果 $X$ 是一个随机变量，它总体分布的均值为 $\mu$ 和标准差(SD)为 $\sigma$ ，而$a$是一个常数，那么 $aX$ 这个新的随机变量分布的均值是$a\mu$，标准差为$|a|\sigma$。为了说明这一点，我们就假设，每只小鼠的体重都减去10，那么平均体重也应该会减去10，类似的，如果我们将小鼠的体重单位由克(g, gram)换成毫克(mg, milligram)，也就相当于把原来的数字都乘以了1000，那么数字的扩散程度（这个扩散程度就是方差或标准差）也会变大。

这就说明了，如果我们获得了 $N$ 个样本，那么这个公式 $\frac{\overline{Y}-\mu}{\sigma_{Y} / \sqrt{N}}$ 就接近于均值为0，标准差为1的正态分布。

现在我们感兴趣的是两个样本的均值的差值。这个数学公式就比较有用了。

如果我们有两个随机变量，即 $X$ 和 $Y$ ，它们的总体均值与总体均值差为$\mu_{X}$和$\mu_{Y}$，$\sigma_{X}$和$\sigma_{Y}$，然后我们就会得到这些结果：

$X+Y$的均值：$\mu_{X}+\mu_{Y}$。

也会得到$Y-X$的均值，其实这也相当于$Y-X=Y+aX$，其中$a=-1$，这其实也就是说$Y-X$的均值为$\mu_{Y}-\mu_{X}$。

如果说$X$和$Y$是相互独立的样本，那么$Y+X$的方差，就是$\sigma_{Y}^2+\sigma_{X}^2$，从这里我们可以推到，$Y-X$，也即$Y+aX$（其中a为-1）的方差是$\sigma_{Y}^2+a^2\sigma_{X}^2=\sigma_{Y}^2+\sigma_{X}^2$。因此说，$Y-X$这个差值的方差也是两个样本的方差和，这有点反直觉，我们需要知道的就是，X与Y这两个样本是相互独立的，它们确实不相互影响。

从数学角度来理解上面的内容对我们的研究目（两个样本的均值，以及均值的差值）。因此这两个样本都服从正态分布，它们的差值也服从正态分布，方差是两个样本总体的方差之和。在零假设（也就是说这两个样本所代表的两个总体的均值没有差异）成立的前提下，样本均值$\bar{Y}-\bar{X}$的差值的分布接近于均值为0（也就是说没有差异），标准差为$\sqrt{\sigma_{X}^{2}+\sigma_{Y}^{2}} / \sqrt{N}$的正态分布，也就是说下面的这个公式：

的值的分布接近于标准正态分布（均值为0，标准差为1）。

通过这种近似计算，我们就能够将计算p值的过程简单化，因为我们已经知道了标准正态分布下任意值范围的比例。例如，这些值中只有5%的比例会大于2（绝对值），计算过程如下所示：

1
pnorm(-2) + (1 - pnorm(2))

计算结果如下所示：

1
2
> pnorm(-2) + (1 - pnorm(2))
[1] 0.04550026

从上面结果我们就可以知道，只有5%的可能性会大于2，我们前面计算的obstdiff的值是3.020833，因此12只小鼠足够了，不需要买更多的小鼠。

但是，这还没完。因为我们并不知道总体的标准差，也就是不知道$\sigma_{X}$与$\sigma_{Y}$。它们是未知的总体参数，但是，我们可以通过样本的标准差对它们进行估算，估算的总体标准差我们不能希腊字母$\mu$来表示，只能用$s$来表示，那么它们的估计值如下所示：

我们需要注意的是，上面的公式里面除以的是$M-1$与$N-1$，而不是前面的$M$与$N$，这有数学方面计算的原因，这里不表，记住就行。但是，为了更加直观地说明一些东西，我们还是要涉及一点。试想，如果你有2个数，把它们画在数轴上，它们之间的距离如果是L，那么中间的那一点就是它们的平均值，这个平均值离它们的距离分别都是0.5L。因此，从这个角度来说，你从这2个数中的一个就能够获取这个信息（这涉及自由度的问题），而不需要2个数，这点不怎么重要，重要的一点是，我们是用了$s_{X}$和$s_{Y}$来估计总体的均值，即$\sigma_{X}$和$\sigma_{Y}$，因为总体均值未知，只能靠样本均值来估计。

因此我们就把前面的那个公式：

定义为下面的这个样子：

如果$M=N$，那么公式就可以化简为如下所示：

CLT告诉我们，当M与N足够大的时候，上面的这个随机变量就会服从均值为0，标准差为1的正态分布。

t分布

CLT的计算依赖于大量的样本，也就是我们指的渐进性结果。当CLT不适用时，还需要另外一种不依赖于渐进性结果的选择。当原始总体来源于一个变量，也就是$Y$时，如果这个$Y$是服从均值为0的标准正态分布，那么我们就可以计算下面这个变量的分布：

这是两个随机变量的比值，因此这个比值不一定服从正态分布。事实上，这个公式的分母会偶然变小，因此就相应地会增加观察到更大值的概率。Gosset最初发现的这种分布，由于他开始在论文中以Student的名义投的稿，因此就称这种分布为Student‘s t-distribution，也就是t分布。

现在我们以小鼠的表型数据为例说明一下，我们创建2向量，一个当作对照总体，一个当作高脂总体，如下所示：

1
2
3
4
5
6
7
8
9
10
url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extd\
ata/mice_pheno.csv"
filename <- "mice_pheno.csv"
download(url,destfile=filename)
dat <- read.csv(filename)
head(dat)
library(dplyr)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist

这里需要注意提，我们假设的是的$y_{1},y_{2},\dots,y_{n}$的分布，而不是随机变量$\bar{Y}$的分布。现在我们来看一下这两个总体的分布：

1
2
3
4
library(rafalib)
mypar(1,2)
hist(hfPopulation)
hist(controlPopulation)

我们可以使用qq图来确认一下上面的分布是比较接近正态分布的，这里我们只作初步了解，更详细的了解在后面。如果在qq图上，点大致都藻在一条直线附近，它们就属于正态分布，如下所示：

1
2
3
4
5
mypar(1,2)
qqnorm(hfPopulation)
qqline(hfPopulation)
qqnorm(controlPopulation)
qqline(controlPopulation)

练习

P47

中心极限定理的实际运用

现在我们使用上面的数据来看一下，中心极限定理是如何逼近样本的均值的：

1
2
3
4
5
6
7
8
9
10
11
12
13
url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extd\
ata/mice_pheno.csv"
filename <- "mice_pheno.csv"
download(url,destfile=filename)
dat <- read.csv(filename)
head(dat)
library(dplyr)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist
mu_hf <- mean(hfPopulation)
mu_control <- mean(controlPopulation)
mu_hf - mu_control

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
> head(dat)
  Sex Diet Bodyweight
1   F   hf      31.94
2   F   hf      32.48
3   F   hf      22.82
4   F   hf      19.92
5   F   hf      32.22
6   F   hf      27.50
> 
> library(dplyr)
> controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
> hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist
> mu_hf <- mean(hfPopulation)
> mu_control <- mean(controlPopulation)
> mu_hf - mu_control
[1] 2.375517

现在计算一下总体的标准差，这里不要使用R的sd()函数，因此这个函数是计算样本的标准差的，它会除以(n-1)，当我们想要对总体进行估计时，按以下代码计算：

1
2
3
4
5
6
7
8
x <- controlPopulation
N <- length(x)
populationvar <- mean((x-mean(x))^2)
populationvar
var(x)
var(x)*(N-1)/N
identical(var(x)*(N-1)/N,populationvar)
identical(var(x), populationvar)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
> populationvar
[1] 11.67205
> var(x)
[1] 11.72416
> var(x)*(N-1)/N
[1] 11.67205
> identical(var(x)*(N-1)/N,populationvar)
[1] TRUE
> identical(var(x), populationvar)
[1] FALSE

因此，为了在数学计算上是准确的，我们不使用sd()或var()函数，我们可以使用rafalib包中的popvar()与popsd()函数，如下所示：

1
2
3
library(rafalib)
sd_hf <- popsd(hfPopulation)
sd_control <- popsd(controlPopulation)

这里需要注意，在实际统计中，我们并不会计算总体参数，我们通常是由样本来估计总体参数，现在我们从总体中抽样，如下所示：

1
2
3
N <- 12
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation,N)

CLT告诉我们，对于更大的N来说，这些样本都近拉于正态分布。作为一种经验，这个N通常是大于30，这仅是一个经验法则，具体来说，它还需要取决于总体的分布。这里我们实际上可以查看一下近似值，并对各种N值进行计算，代码如下所示：

1
2
3
4
5
Ns <- c(3, 12, 25, 50)
B <- 10000
res <- sapply(Ns, function(n){
  replicate(B, mean(sample(hfPopulation, n)) - mean(sample(controlPopulation, n)))
})

现在使用qq图看一下CLT理论对这些数据的近似值的适应情况。如果这些数据非常接近正态分布，那么这些数据点就会落在一个直线上（这个直接是正态分布的分位数(quantiles))。越偏离，说明数据越不符合正态分布，如下所示：

1
2
3
4
5
6
7
8
9
mypar(2,2)
for (i in seq(along=Ns)){
  titleavg <- signif(mean(res[,i]),3)
  titlesd <- signif(popsd(res[,i]),3)
  title <- paste0("N=", Ns[i], "Avg=", titleavg, " SD=", titlesd)
  qqnorm(res[,i], main = title)
  qqline(res[,i], col=2)
  
}

从上面结果可以看出来，第3组数据的拟合结果最合适，这是因为它的总体相对最接近正态分布，平均值也最接近于正态分布。在实际中，我们可以计算这样一个比值：除以估计的标准误差来判断，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Ns <- c(3, 12, 25, 50)
B <- 10000
computetstat <- function(n){
  y <- sample(hfPopulation, n)
  x <- sample(controlPopulation, n)
  (mean(y)- mean(x))/sqrt(var(y)/n+var(x)/n)
}
res <- sapply(Ns, function(n){
  replicate(B, computetstat(n))
})
mypar(2,2)
for (i in seq(along=Ns)){
  qqnorm(res[,i], main = Ns[i])
  qqline(res[,i], col=2)
  
}

从上面可以看出来，当$N=3$时，CLT并不会提供一个有用的估计，当$N=12$时，在较大值的点方面，稍微有点偏离。当$N=25$或$N=50$时，所有的点基本上都在直线上。

也就是说，在这个案例中，只要$N=12$就能证明CLT，就像前言提到的那样，在多数情况下，这种模拟并不会很好。我们在这里只是通过这种模拟手段来说明CLT背后的思想与局限。

练习

P52

t检验的实际计算

现在我们来看一下如何通过计算得到p值，先来导入数据，在这一步中，我们要确定哪些是treatment组，哪些是control组，并且计算出它们的均值差异，如下所示：

1
2
3
4
5
6
7
8
library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)
control <- filter(dat,Diet=="chow") %>% select(Bodyweight) %>% unlist
treatment <- filter(dat,Diet=="hf") %>% select(Bodyweight) %>% unlist
diff <- mean(treatment) - mean(control)
print(diff)

计算结果如下所示：

1
2
> print(diff)
[1] 3.020833

现在我们开始计算p值，上面的代码中有一个diff变量，可以称为observed effect size，这是一个随机变量，我们还知道，在零假设的前提下，diff均值的分布是0，而这个随机变量分布的标准误则是总体的标准差(population standard deviation)除了样本数目的平方根，如下所示：

我们使用样本样本标准差(sample standard deviation)来对总体的标准差(population standard deviation)进行估计，在R中，也可以使用的sd()函数来计算标准误(SE)，如下所示：

1
sd(control)/sqrt(length(control))

结果如下所示：

1
2
> sd(control)/sqrt(length(control))
[1] 0.8725323

这个值就是样本均值(sample average)的SE，不过实际上，我们相想要的是diff的SE。我们前面已经知道了，统计学理论告诉我们，两个随机变量差的方差是原两个随机变量方差的和，因此我们计算一下方差以及平方根：

1
2
3
4
5
se <- sqrt(
var(treatment)/length(treatment) +
var(control)/length(control)
)
se

结果如下所示：

1
2
> se
[1] 1.469867

统计学理论还告诉我们，如果我们将一个随机变量除以其SE，我们就会得到一个SE为1的新的随机变量，如下所示：

1
2
tstat <- diff/se
tstat

结果如下所示：

1
2
3
> tstat <- diff/se
> tstat
[1] 2.055174

上面计算出来的tstat就是我们称的t统计量(t-statistic)。这是两个随机变量的比值，因此它也是一个随机变量。一旦我们知道了这个随机的分布，我们就很容易计算出其p值。

前面我们提到，CLT可以告诉我们，如果有大样本(large sample size)，那么样本的均值mean(treatment)和mean(control)都服从正态分布。统计学理论告诉我们，两个正态分布的随机变量的差值还服从正态分布，因此CLT也告诉我们，tstat接近于均值为0(零假设)，SD为1（我们除以它的SE）的正态分布。

现在我们要计算一下p值，此时我们需要问的一个问题就是：一个服从正态分布的随机变量有多大的概率会大于diff？R语言中的有一个pnorm()函数有解决这个问题。pnorm(a)会以返回一个值，这个值就是概率，它表示在一个标准正态分布中，一个随机变量低于a 的概率，具体这个函数的原理与用法，可以参考这篇笔记《R语言中dnorm, pnorm, qnorm与rnorm以及随机数》。

如果要计算大于a的概率，那么可以使用1-pnorm(a)，如果我们想知道看到一些像diff这种极端事件的概率：例如要么小于（更多的时候指的是负值）-abs(diff)，要么是大于abs(diff)，那么我们也可以过计算这两个尾部区域的概率：

1
2
3
4
righttail <- 1-pnorm(abs(tstat))
lefttail <- pnorm(-abs(tstat))
pval <- lefttail + righttail
print(pval)

计算结果如下所示：

1
2
> print(pval)
[1] 0.0398622

在这个案例中，p值小于0.05，根据传统的p值阈值0.05，我们就可以下给结论，这种差值有统计学意义。

现在我们就面临一个问题，CLT只有在大样本的前提下才有效，12个样本这个数目够么？根据经验法则，通常大于30个样本（这仅仅是经验法则），就能满足CLT。我们计算的这个p值，只有在假设成立的前提下才是有效的近似值，而实际情况并非如此，还好，除了CLT之外，我们还有另外一种选择。

t分布的实际计算

如果一个总体的分布是正态分布，那么我们在不需要CLT思想的前提下，就可以计算出t-statistic的精确分布。但是，如果我们只有小样本，我们就很难知道总体是否是正态的。但是，对于一些常见的统计量，例如体重，我们就可以猜测它的总体分布非常接近于正态分布，因此我们可以使用这种近似。此外，我们可以使用对样本画qq图，它会显示出样本的分布，如下所示：

1
2
3
4
5
6
library(rafalib)
mypar(1,2)
qqnorm(treatment)
qqline(treatment,col=2)
qqnorm(control)
qqline(control,col=2)

如果我们使用这种近似计算，那么从统计学理论角度来讲我们随机变量tstat的分布就服从一个t分布(-t-distribution)，这种分布比正态分布更复杂。t分布没有像正态分布那样的位置参数( location parameter )，正态分布的参数是自由度(degrees of freedom)，R中的t.test()函数可以用于计算相应的p值，如下所示：

1
t.test(treatment, control)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
> t.test(treatment, control)
	Welch Two Sample t-test
data:  treatment and control
t = 2.0552, df = 20.236, p-value = 0.053
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.04296563  6.08463229
sample estimates:
mean of x mean of y 
 26.83417  23.81333

如果只想看p值，那么使用以下代码即可，如下所示：

1
2
result <- t.test(treatment,control)
result$p.value

结果如下：

1
2
3
> result <- t.test(treatment,control)
> result$p.value
[1] 0.05299888

从计算结果来看，这个p值有点大，这是因为在CLT近似的计算中，tstat的分布实际上是固定的（对于大样本来说，分母实际上就是固定的），而t分布的近似计算则是考虑到了分母（差值的标准误）是一个随机变量，是不固定的，样本数目越小，那分母的变化就越大。

这样来看的，我们使用了两种方法，得到了2个不同的p值，这在数据分析过程中很常见，因为我们是使用了不同的前提，不同的近似，因此会得到不同的结果。在后面讲到的功效计算(power calculation)中，我们会讲到I类错误，II类错误。在这里我们只说一下，使用CLT近似会错误地拒绝( incorrectly reject )零假设(假阳性)，而t分布则会错误地接受(inccorrectly acept)零假设(假阴性)。

在R中计算t检验如下所示：

1
2
3
4
5
6
7
library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
control <- filter(dat,Diet=="chow") %>% select(Bodyweight)
treatment <- filter(dat,Diet=="hf") %>% select(Bodyweight)
t.test(treatment,control)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
> t.test(treatment,control)
	Welch Two Sample t-test
data:  treatment and control
t = 7.1932, df = 735.02, p-value = 1.563e-12
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 2.231533 3.906857
sample estimates:
mean of x mean of y 
 30.48201  27.41281

置信区间(Confidence Intervals)

在生命科学研究中，仅在研究结论中报道p值是不够的，因为统计学上的显著性无法保证科学研究上的意义。例如在对体重的差异进行统计时，使用大样本的情况下，即使是1毫克的差异，也有可能有统计学意义，不过这种1毫克的差异是一个重要的发现么？我们能否说明，某个因素导致了不到1%的体重变化？仅报道p值是无法提供一些有价值的信息的，也不会提供有关效应量(the effect size)的信息，效应量(effect size)就是我们前面提到的观察到的差异(observed difference)。有的时候，效应量会用除以对照组的均值，因此会表示为百分比。

为了纠正仅报道p值偏差，另外一个统计量就提了出来，即置信区间(confidence intervals)。置信区间包括了估计的效应量(estimated effect size)，以及与这个估计有关的不确定性(uncertainty)，现在我们使用这批小鼠数据来计算一下置信区间。

总体均值的置信区间

在我们展示如何过计算两组差值的置信区间之前，我们先来看一下如何计算对照组雌性小鼠的总体均值的置信区间。随后，我们会计算样本的置留区间，如下所示：

1
2
3
4
5
6
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
chowPopulation <- dat[dat$Sex=="F" & dat$Diet=="chow",3]
mu_chow <- mean(chowPopulation)
print(mu_chow)

计算结果如下所示：

1
2
> print(mu_chow)
[1] 23.89338

现在计算的结果就是 总体均值$\mu_{X}$。

现在我们对这个结果进行估计，在实际分析过程中，我们并不会获取所有的总体，因此我们也像前面计算p值那样，我们来演示一下，通过抽样的方法来计算置信区间，现在我们先抽30个样本，如下所示：

1
2
3
N <- 30
chow <- sample(chowPopulation, N)
print(mean(chow))

计算结果如下所示：

1
2
> print(mean(chow))
[1] 23.95033

我们知道这个均值是一个随机变量，因此样本均值并不是一个非常好的估计。实际上，我们使用这个案例只是为了计算这个样本的均值，每次计算并非完全一样，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> N <- 30
> chow <- sample(chowPopulation, N)
> print(mean(chow))
[1] 23.48467
> 
> N <- 30
> chow <- sample(chowPopulation, N)
> print(mean(chow))
[1] 25.22967
> 
> N <- 30
> chow <- sample(chowPopulation, N)
> print(mean(chow))
[1] 23.732

如上所示，每次抽样计算的均值都不一样。现在我们来了解一下置信区间，置信区间一种报道你数据中样本均值的一种统计学表示方法，它能明确地显示出随机变量的变异程度。

现在我们的的样本是30，CLT告诉我们，$\bar{X}$或chow组的均值(mean)服从一个均值为$\mu_{X}$（这个是chowPopulaiton这个总体的均值，即mean(chowPopulation)，这个均值服从均值为$\mu_{X}$，标准误$s_{X} / \sqrt{N}$的标准分布，计算过程如下所示：

1
2
se <- sd(chow)/sqrt(N)
print(se)

计算结果如下所示：

1
2
3
> se <- sd(chow)/sqrt(N)
> print(se)
[1] 0.597656

定义区间

95%置信区间是一个随机区间，这个区间有95%概率落在我们估计值的参数中。这里需要注意的是，说95%的随机区间将会包括真值，并不等于说真值有95%的概率落在我们的区间中。为了构建置信区间，我们需要注意，CLT告诉我们，$\sqrt{N}\left(\overline{X}-\mu_{X}\right) / s_{X}$服从均值为0，SD为1的正态分布，这个概率也就是以下事件的概率：

在R中的计算为：

1
pnorm(2) - pnorm(-2)

计算结果如下所示：

1
2
> pnorm(2) - pnorm(-2)
[1] 0.9544997

这个值与95%很接近，也就是qnorm(1-0.05/2)的这个值，现在我们把上面的公式变换一下，把$\mu_{X}$放中间，就成了如下的样子：

也就是上面这个值的概率是95%。需要注意的是，区间$\overline{X} \pm 2 s_{X} / \sqrt{N}$是两个边界，它们是随机的。此外，我们需要明确一下，置信区间的定义是，95%的随机区间（random intervals)包含真正的固定值$\mu_{X}$。对于一个计算结果来说，计算出来的这个特定的区间，它要么包括固定的总体均值$\mu$，要么不包括。

现在我们来模拟一下其中的思路，如下所示：

1
2
3
4
5
6
7
N <- 30
chow <- sample(chowPopulation, N)
se <- sd(chow)/sqrt(N)
Q <- qnorm(1 - 0.05/2)
interval <- c(mean(chow)-Q*se, mean(chow) + Q*se)
interval
interval[1] < mu_chow & interval[2] > mu_chow

结果如下所示：

1
2
3
4
> interval
[1] 22.17089 24.58044
> interval[1] < mu_chow & interval[2] > mu_chow
[1] TRUE

从上面我们可以看到，这个区间覆盖了总体均值$\mu_{X}$或mean(chowPopulation)。但是，如果我们再取另外一批样本的话，有可能不会覆盖总体均值$\mu_{X}$，但是，统计学理论告诉我们，经过均数次抽样后，这个样本覆盖均值的概率是95%。由于我们能够获取总体数据，现在我们来看一些样本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(rafalib)
B <- 250
mypar()
plot(mean(chowPopulation)+c(-7, 7), c(1,1), type="n",
     xlab="weight", ylab = "interval", ylim=c(1,B))
abline(v=mean(chowPopulation))
for (i in 1:B){
  chow <- sample(chowPopulation, N)
  sd <- sd(chow)/sqrt(N)
  interval <- c(mean(chow)-Q*se, mean(chow)+Q*se)
  covered <- 
    mean(chowPopulation) <= interval[2] & mean(chowPopulation) >= interval[1]
  color <- ifelse(covered, 1,2)
  lines(interval, c(i, i), col= color)
}

在这个案例中，我们模拟了250次随机抽样，计算了250次置信区间，其中颜色表示这个区间是否包含总体均值，其中绿色是包括的，棕色是不包括的。

上面这个案例的置信区间比较大（这个案例除以的是是$\sqrt{5}$，而非$\sqrt{30}$），我们可以看到很多区间并未覆盖$\mu_{X}$，这是因为CLT错误地告诉了我们(chow)的均值是近似地接近于正态分布，实际上这个数据是不太符合正态分布的（在接近$\pm \infty$的地方，它有一个比较扁平的尾巴）。这种错误会影响我们计算Q值，而这个Q值在计算的时，是以正态分布为前提,，使用qnorm()来计算的。在这种情况下，数据会更符合t分布，现在我们使用qt()函数来计算Q值，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
mypar()
plot(mean(chowPopulation)+ c(-7,7), c(1,1),type="n",
     xlab="weight",ylab="interval",ylim=c(1,B))
abline(v=mean(chowPopulation))
Q <- qt(1-0.05/2, df=4)
N <- 5
for (i in 1:B){
  chow <- sample(chowPopulation, N)
  se <- sd(chow)/sqrt(N)
  interval <- c(mean(chow)-Q*se, mean(chow)+Q*se)
  covered <- mean(chowPopulation) <= interval[2] & mean(chowPopulation) >= interval[1]
  color <- ifelse(covered, 1,2)
  lines(interval, c(i,i), col=color)
}

上面这个图中，我们使用了小样本来模拟了250次计算，得到了它的95%置信区间，这个置信区间是基于t分布进行计算的。与前面的那个案例相比，这个案例中置信区间更大，这是因为t分布的尾部更为扁平，得到的置信区间就更大，现在我们来计算一下t分布的累积分分布，并拿它与正态分布比较一下，如下所示：

1
2
qt(1-0.05/2, df=4)
qnorm(1-0.05/2)

计算结果如下所示：

1
2
3
4
> qt(1-0.05/2, df=4)
[1] 2.776445
> qnorm(1-0.05/2)
[1] 1.959964

从上面结果就可以看出来，t分布的置信区间更大，因此它才能以更大的概率(例如95%)来覆盖$\mu_{X}$。

置信区间与p值的关系

作者推荐报道实验结果时使用置信区间，而非仅用p值，对于某些原因，我们可能只需要提供p值即可，或者是提供0.05或0.01水平上的统计结果，其实置信区间也能提供这些信息。

当我们谈论t检验(t-test)的p值时，我们实际上是在谈论我们观察到的这两个样本的差值($\overline{Y}-\overline{X}$)以及更极端情况下差值的概率，就像是谈论两个总体之间的总体均值之差等于0这种极端情况下的概率。因此我们可以构建一个含有这个观察到的差值的置信区间，不过我们不再写为$\overline{Y}-\overline{X}$这种形成，而是写为一种新的形式，即$d \equiv \overline{Y}-\overline{X}$。

当我们来写一个差值的95%置信区间时，使用CLT则会写为$d \pm 2 s_{d} / \sqrt{N}$这种形式，这个结果并不包含0（假阳性）。因为置信区间不含0，这就暗示了，$d-2 s_{d} / \sqrt{N}>0$或$d+2 s_{d} / \sqrt{N}<0$，这种表达形式还说明，$\sqrt{n} d="" s_{d}="">2$或$\sqrt{N} d / s_{d}<2$。这说明，t统计量比2更极端，反过来又说明，p值必然小于0.05（为了更加精确地计算，可以使用qnorm(0.05/2)来替换2）。如果使用t分布来取代CLT，则使用使用qt(0.05/2, df=N-2。总之，一个95%或99%置信区间并不包含0，p它们的p值必然小于0.05或小于0.01。

现在我们使用t检验来计算一下d的置信区间，如下所示：

1
t.test(treatment,control,conf.level=0.95)$conf.int

结果如下所示：

1
2
3
4
> t.test(treatment,control,conf.level=0.95)$conf.int
[1] 2.231533 3.906857
attr(,"conf.level")
[1] 0.95

如果我们把置信水平改为0.9，如下所示：

1
2
3
4
> t.test(treatment,control,conf.level=0.90)$conf.int
[1] 2.366479 3.771911
attr(,"conf.level")
[1] 0.9

可以发现，置信区间变小了。

功效计算

在前面的案例中，我们研究了不同包含对小鼠体重的影响，由于我们在使用这个案例的时候，我们研究的是总体，发现了这两种包含对小鼠的体重有着明显的影响，计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% dplyr::select(Bodyweight) %>% unlist
mu_hf <- mean(hfPopulation)
mu_control <- mean(controlPopulation)
print(mu_hf - mu_control)
print((mu_hf - mu_control)/mu_control*100)

计算结果如下所示：

1
2
3
4
> print(mu_hf - mu_control)
[1] 2.375517
> print((mu_hf - mu_control)/mu_control*100)
[1] 9.942157

某些情况下，我们会采用抽样的形式，抽取一定数目的小鼠体重作为样本，此时就不能使用总体的计算方法，而是采用t检验，此时计算得到的p值不一定小于0.05，例如，我们随机抽取5只小鼠来计算一下，如下所示：

1
2
3
4
5
set.seed(1)
N <- 5
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation, N)
t.test(hf,control)$p.value

计算结果如下所示：

1
2
> t.test(hf,control)$p.value
[1] 0.1410204

前后两种的计算方法出现了差异，第一次计算（使用小鼠的总体来进行计算）与第二次计算（随机抽两组的5只小鼠来计算）的p值不一样，是否出现了错误？

在第二次计算中，p值大于0.05，因此我们不拒绝零假设，从而说饮食对小鼠的体重没有影响么？

能否这么讲？答案不是能。

我们只能说，无法拒绝零假设（注意，没有后面那半句：包含对小鼠的体重没有影响）。我们无法拒绝零假设，并不意味着零假设一定为真。在这个特殊的案例中，这里我们面临的一个问题就是，我们没有足够的功效(power)。如果你从事科学研究，很有可能在某些时候不得不做一个功效计算。在很多情况下，这是一种首选义务，因此它可以帮助你在人的研究中，避免使用过多的小鼠，或者是避免使用过多的受试患者。现在我们就来解释一下统计功效(statistical power)的含义。

错误类型

我们要有这么一个思路，一旦我们做统计，就有可能出现计算错误，这也就是为什么p值不等于0。在进行统计时，当p值很小时，这说明我们可以拒绝零假设，但是这种拒绝也有可能犯错误（虽然犯错误的概率很低），例如零假设确实是真时。例如当p值等于0.05时，那么我们重复20实验，就有可能1次会发生零假设成功的事件，统计学家把这种错误称为I类错误(type I error)。

I类错误的定义就是，零假设成立时，我们拒绝了（其实不应该拒绝的），也就是所谓假阳性。这里为什么使用0.05，而不使用0.00001呢？这是因为如果阈值设得过低，成本会非常高。此时我们还要引入II类错误(type I error)，II类错误是指，零假设不成立，但是我们接受了零假设，也就是所谓的假阴性。前面我们抽样了5只小鼠，计算得到的值大于0.05，因此我们无法拒绝零假设（在0.05这个水平上），此时所犯的错误就是II类错误，也就是假阴性（从总体上看，实际上是有差异的，但是我们没有发现）。如果我们把阈值提高到0.25，我们就不会犯这个错误了（计算出的p值为0.1410204），不过通常情况下，不会这么干。

0.05与0.01阈值(cut-off)

多数的期刊或监管机构都坚持0.01或0.05的显著性水平，当然这么做无可厚非。我们本书的目标就之一就是让读者对p值以及置信区间有一个更深入的理解，以便读者能够在一些情况下做出合理的选择。不幸的是，这些阈值在某些程度上已经滥用，不过这是又是一个容易引起争论的话量，此处不表。

功效计算

功效(power)是指：当零假设为假时，拒绝它的概率，原书的描述如下所示：

Power is the probability of rejecting the null when the null is false.

不过“零假设为假”这种情况比较复杂，因为在许多情况下，零假设都有可能为假。$\Delta \equiv \overline{Y}-\overline{X}$可以是任意的，功效实际上要依赖于这个参数，功效同时还依赖于你估计值的标准误，反过来，标准误又依赖于样本数目，以及总体标准差。实际情况中，我们并不知道一些$\Delta$，$\sigma_{X}$，$\sigma_{Y}$值以及样本数目这些值。但是通过统计学理论，我们可以计算功效，R中的pwr包可以进行这样的计算，现在我们来模拟一下这个过程：

第一步：确定样本数目，我们确定为12，即N <- 12；

第二步：设定拒绝零假设的阈值，这里为0.05，也就是常见的p值水平，即alpha <- 0.05；

第三步：设计重复抽样次数，我们会重复地进行抽样，来计算一下拒绝零假设的比例，这里设为2000，即B <- 2000；

这个模拟过程就是：我们会从对照组和处理组中抽取N个样本，使用t检验来计算p值，观察p值是否小于0.05，现在我们把这个过程写为一个函数，代码如下所示：

1
2
3
4
5
6
7
8
9
10
N <- 12
alpha <- 0.05
B <- 2000
reject <- function(N, alpha=0.05){
  hf <- sample(hfPopulation, N)
  control <- sample(controlPopulation, N)
  pval <- t.test(hf, control)$p.value
  pval < alpha
}

如果样本数目为12，那么我们是否能拒绝零假设呢，如下所示：

1
2
> reject(12)
[1] FALSE

答案是不能。现在我们使用replicate()函数来重复2000次，如下所示：

1
2
3
N <- 12
rejections <- replicate(B, reject(N))
mean(rejections)

计算结果如下所示：

1
2
> mean(rejections)
[1] 0.218

其中rejecctions的结果是一个逻辑值，即里面只有TRUE或FALSE，其中TRUE为1，FALSE是0，那么我们计算的均值，即mean(rejections)就是TRUE所占的比例，因此我们计算出来的功效为0.218。这就是为什么，当我们知道零假设为假时，t检验也无法拒绝的原因。当样本数目为12时，功效是21.8%，对为能够在0.05水平上降低假阳性，那么我们可以设置更高的阈值，但这会导致II类错误的增多。

现在我们来看一下，功效是如何随着N的变化而变化的，如下所示：

1
2
3
4
5
6
Ns <- seq(5,50, 5)
power <- sapply(Ns, function(N){
  rejections <- replicate(B, reject(N))
  mean(rejections)
})
power

计算结果如下所示：

1
2
> power
 [1] 0.0960 0.1640 0.2775 0.3775 0.4795 0.5540 0.6415 0.7075 0.7725 0.8160

从上面的结果我们可以发现，随着N的增加，功效也在增加，曲线如下所示：

1
plot(Ns, power, type = "b")

再来看一个案例。我们把alpha这个拒绝阈值改变一下，再来看一下功效的变化。如果想要降低I类错误，那么功效就会越小，也就说，我们需要在两类错误之间进行权衡，现在我们看下面的代码，在保证N值相同的前提下，改变alpha后功效的大小：

1
2
3
4
5
6
7
N <- 30
alphas <- c(0.1, 0.05, 0.01, 0.001, 0.001)
power <- sapply(alphas, function(alpha){
  rejections <- replicate(B, reject(N, alpha=alpha))
  mean(rejections)
})
plot(alphas, power, xlab="alpha", type="b", log="x")

变化曲线如下所示：

在实际统计中，没有绝对的功效，也没有绝对的alpha值，但是我们要理解这背后的统计学硬。

替代假设下的p值

当零假设为假，替代假设为真时，那么p值有的时候就比较任意了。当替代假设为真时，我们通过增加样本数目，可以得到一个尽可能小的p值，通过从总体中不断抽取更大的样本数，我们就会知道p值的这种特性，现在看一下这个过程。

第一步：构建一个函数，这个函数能够返回一个样本数目为N时的p值，如下所示：

1
2
3
4
5
calculatePvalue <- function(N){
  hf <- sample(hfPopulation,N)
  control <- sample(controlPopulation, N)
  t.test(hf, control)$p.value
}

第二步：设定样本数目，对于每个样本数目的计算，我们会得到一系列p值，如下所示：

1
2
Ns <- seq(10, 200, by=10)
Ns_rep <- rep(Ns, each=10)

第三步：计算p值，如下所示：

1
pvalues <- sapply(Ns_rep, calculatePvalue)

第四步：绘制样本数目与其对应的p值，如下所示：

1
2
3
plot(Ns_rep, pvalues, log="y", xlab="Sample size",
     ylab="p-value")
abline(h=c(0.01, 0.05), col="red", lwd=2)

上图是随着样本数目的变化，p值的变化。从上面结果可以看出，当替代假设成立时，随着样本数目的增多，p值会越来越小，从上面的两条红线分别是0.01与0.05的区间。

一旦我们在某个阈值上拒绝了零假设，我们如果想要获得一个更小的p值，那么就需要更多的样本，例如实验小鼠。样本增大会增加我们对差值$\Delta$估计的精确性，但是，这个p值的降低仅仅是数学计算的自然结果。因为在计算p值时，它公式中的分母是样本数目的平方根，即$\sqrt{N}$，即如果$\Delta$不等于0，随着N的增加，p值必然降低。

因此，一个好的统计学结果需要列出效应量(effect size)与置信区间。效应量的计算是：将差值(difference)和置信区间(confidence interval)除以对照总体的均值，获得一个百分数，如下所示：

1
2
3
4
5
6
N <- 12
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation, N)
diff <- mean(hf) - mean(control)
diff/mean(control)*100
t.test(hf, control)$conf.int / mean(control) * 100

计算结果为：

1
2
3
4
5
6
> diff/mean(control)*100
[1] 5.783149
> t.test(hf, control)$conf.int / mean(control) * 100
[1] -8.169363 19.735661
attr(,"conf.level")
[1] 0.95

此外，我们还可以写出一个名为Cohen's d的统计量，它是两组之间的差值除以这两组总和的标准差，如下所示：

1
2
sd_pool <- sqrt(((N-1)*var(hf) + (N-1)*var(control))/(2*N - 2))
diff / sd_pool

计算结果如下所示：

1
2
3
> sd_pool <- sqrt(((N-1)*var(hf) + (N-1)*var(control))/(2*N - 2))
> diff / sd_pool
[1] 0.3510624

这个结果告诉我们，hf组的小鼠体重的均值与对照组(control)小鼠体重平均值偏离多少个标准差。在替代假设下，t统计量会随着样本数目的增加而增加，但是效应量(effect size)与Cohen's d却会变得更加精确。

练习

P76